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The phase behavior of charge-stabilized colloidal suspensions is modeled by a combination of 
response theory for electrostatic interparticle interactions and variational theory for free energies. 
Integrating out degrees of freedom of the microions (counterions, salt ions), the macroion-microion 
mixture is mapped onto a one-component system governed by effective macroion interactions. Linear 
response of microions to the electrostatic potential of the macroions results in a screened-Coulomb 
(Yukawa) effective pair potential and a one-body volume energy, while nonlinear response modifies 
the effective interactions [A. R. Denton, Phys. Rev. E TO, 031404 (2004)]. The volume energy and 
effective pair potential are taken as input to a variational free energy, based on thermodynamic 
perturbation theory. For both linear and first-order nonlinear effective interactions, a coexistence 
analysis applied to aqueous suspensions of highly charged macroions and monovalent microions yields 
bulk separation of macroion-rich and macroion-poor phases below a critical salt concentration, in 
qualitative agreement with predictions of related linearized theories [R. van Roij, M. Dijkstra, and 
J.-P. Hansen, Phys. Rev. E 59, 2010 (1999); P. B. Warren, J. Chera. Phys. 112, 4683 (2000)]. It 
is concluded that nonlinear screening can modify phase behavior but does not necessarily suppress 
bulk phase separation of deionized suspensions. 



I. INTRODUCTION 

Mounting evidence from a variety of experiments suggests that colloidal suspensions 0, 0> of highly 
charged macroions and monovalent microions (counterions and coions) can separate into macroion-rich and 
-poor bulk phases at low salt concentrations. Reported observations - in aqueous suspensions at sub- 
millimolar ionic strengths - describe liquid- vapor coexistence , stable voids [H IS 13 j contracted crystal 
lattices H, H, 0] , and metastable crystallites 0] . Such phenomena suggest an unusual form of interparticle 
cohesion, inconsistent with the long-ranged repulsive electrostatic pair interactions that prevail at low ionic 
strengths |l 2l. a nd in apparent conflict with the classic theory of Derjaguin, Landau, Verwey, and Overbeek 
(DLVO) [13, Hj]) which so successfully describes phase stability with respect to coagulation at higher salt 
concentrations. Observations of bulk phase separation in deionized suspensions are therefore often considered 
anomalous. 

Reports of anomalous phase behavior in charged colloids have been variously disputed |l5j |. at- 
tributed to impurities JlBj. |17[, or interpreted as genuine manifestations of like-charge interparticle at- 
traction 0, H, S S IE 1st lid llfl. whether pairwise or many-body in origin 0, 0, |2(j ■ Although some 
particle-tracking experiments [ill l2lL I22I |23| appear to exhibit attractive forces between isolated pairs of 
tightly confined macroions, recent studies, based on refined optical imaging methods, have found no attrac- 
tion [24(. Furthermore, mathematical proofs that Poisson-Boltzmann theory predicts purely repulsive pair 
interactions |2a l26ll27| relegate any possible pair attraction to the influence of counterion correlations, ne- 
glected by the mean-field theory. It is now widely accepted that correlations among multivalent counterions 
can induce attraction between like-charged surfac es I2& 129.1301 l3ll l32| , as well as condensation of DNA and 
other polyelectrolytes HHH HI |M EUS ImEoTlTlV The key issue motivating the present study 
is whether relatively weakly correlated monovalent counterions can similarly destabilize deionized colloidal 
suspensions. 

Further evidence for effective attractive interactions in charged colloids comes from computer simulations. 
Monte Carlo simulations 0> HE El C3 of the primitive model of asymmetric electrolytes - macroions 
and microions, in a dielectic continuum, directly interacting via repulsive Coulomb pair potentials - ex- 
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hibit macroion attraction and instabilities toward macroion aggregation at high electrostatic couplings. 
Short-ranged attractions have been linked to s pati al correlations among counterions localized near different 
macroions 0, E3 , or to Coulomb depletion |44j , while long-ranged attractions have been attributed to 
overcharging of macroions [45| . System parameters thus far explored correspond to relatively strongly cor- 
related (multivalent) counterions and relatively small macroion-to-counterion size and charge asymmetries. 
Computational advances, however, are rapidly closing the gap that currently prevents direct comparison of 
simulations and experiments. 

Many theoretical studies of interparticle interactions and phase behavior in charged colloids have been mo- 
tivated by the puzzling results of experiments and simulations. Among various analytical and computational 
approaches, recently reviewed Eg, Ell UM , are integral-equation, Poisson-Boltzmann, density-functional, 
Debye-Huckel, and response theories. In seminal work, van Roij et al. El H2 , described the phase be- 

havior of charged colloids within an effective one-component model governed by density-dependent effective 
interactions. Combining a linearized density- functional theory |53| for the effective pair and one-body (vol- 
ume energy) potentials with a variational theory for the free energy, these authors predicted counterion-driven 
bulk phase separation in deionized suspensions of highly charged macroions below a critical salt concentra- 
tion. Subsequently, Warren |54j applied an extended Debye-Huckel (linearized Poisson-Boltzmann) the ory 
and predicted similarly unusual phase separation at low salt concentrations. Statistical mechanical |55l |56( 
and linear-response |57l l58l l59| methods, based on closely related linearization approximations, yield similar 
effective electrostatic interactions. 

Several recent studies, based on Poisson-Boltzmann cell models [(KjEniE^] an d extensions of Debye-Hiickel 
theory [6i| , have suggested that predicted instabilities of charged colloids towards phase separation may be 
mere artifacts of linearization. The main purpose of the present study is to directly test this suggestion by 
explicitly calculating the effect of nonlinear screening on the phase behavior of charged colloids. Working 
within the framework of the effective one-component model and response theory |57l IHsl IHflj . we input 
nonlinear corrections to the effective pair potential and volume energy into an accurate variational free 
energy and analyze thermodynamic phase behavior. The central conclusion of the paper is that nonlinear 
effects can modify phase behavior of deionized suspensions, but do not necessarily suppress counterion-driven 
phase separation. 

Outlining the remainder of the paper, Sec. [H] first defines the model colloidal suspension. Section HTD next 
reviews the response theory for effective interactions and describes a variational perturbation theory for 
the free energy. Section IIVI presents and discusses numerical results - most importantly, equilibrium phase 
diagrams obtained from a coexistence analysis. Finally, Sec. IVl summarizes and concludes. 



II. MODEL SYSTEM 



The system of interest comprises N m negatively charged colloidal macroions, N c positively charged coun- 
terions, and N s pairs of oppositely charged salt ions all dispersed in a solvent. The macroions are modeled as 
charged hard spheres of radius a (diameter a) and effective valence Z , as depicted in Fig. The macroion 
surface charge — Ze is best interpreted as an effective (renormalized) charge, equal to the bare charge less 
the combined charge of any strongly associated counterions. The effective charge is assumed to be uniformly 
distributed over the surface and fixed, independent of thermodynamic state. The counterions and salt ions 
are modeled as point charges of valence z, whose number N c is determined by the condition of overall charge 
neutrality: ZN m = zN c . Numerical results are presented below fSec. HVII for the case of monovalent (z = 1) 
microions. The microions number N + = N c + N s positive and N- — N s negative, totaling — N c + 2N S . 

Working within the primitive model of charged colloids, we approximate the solvent as a dielectric con- 
tinuum, characterized entirely by a dielectric constant e. We further assume a rigid-ion model, ignoring 
van der Waals 12] and polarization l65l l6^| interactions, which are dominated by longer-ranged direct 
electrostatic interactions at low ionic strengths. The system is imagined to be in thermal equilibrium with 
a heat bath at constant temperature and in chemical (Donnan) equilibrium with a salt reservoir (e.g., via a 
semi-permeable membrane or ion-exchange resin), which fixes the salt chemical potential. Having specified 
the model system, we turn next to methods for describing electrostatic interactions and thermodynamic 
phase behavior. 
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III. METHODS 
A. Effective Electrostatic Interactions 

1. One-Component Mapping 

Response theory of effective interactions is fundamentally based on mapping a multi-component mixture 
onto a one-component system governed by an effective Hamiltonian [671 . When applied to charged colloids, 
polyclcctrolyt.es. and other ionic systems, the mapping involves integrating out from the partition function 
the degrees of freedom of the microions ■ The resulting effective interactions between macroions depend 
on the perturbation of the microion distribution by the "external" potential of the macroions. The response 
of the microions to the macroions is linear |57ll58j for suspensions of weakly charged macroions, but becomes 
increasingly nonlinear |59j | as the macroion valence increases and as the salt concentration decreases. Here 
we briefly review the theory, referring the reader to refs. [57|, |5^, IHjJ for further details. 

In the simplest case of a salt-free suspension, the Hamiltonian may be expressed as 

H = H m ({R}) + H c ({r}) + H mc ({R}, {r}), (1) 
where {R} and {r} denote coordinates of macroions and microions, respectively, 

H m = H m ({R}) + - « mm (|R<-R,-|) (2) 

is the Hamiltonian of the macroions alone, 7?hs is the Hamiltonian of a hard-sphere (HS) system, v mm (r) = 
Z 2 e 2 /er, r > er, is the bare Coulomb pair potential between macroions, 



1 N ° 

Hc = K c + - ]T Mki-'il) (3) 

is the counterion Hamiltonian, K c is the counterion kinetic energy, v cc {r) = z 2 e 2 /er is the pair potential 
between counterions, 

N m N c 

H mc = ^2^2v mc (\R i -r j \) (4) 

i=l j=l 

is the total macroion-counterion interaction energy, and v mc (r) — Zze 2 / er, r > a, is the macroion-counterion 
pair potential. 

The mapping from the macroion-counterion mixture to an effective one-component system of pseudo- 
macroions begins with the canonical partition function 

Z = «exp(-/3ff)) c ) m , (5) 

where f3 = l/(fcgT) at temperature T and angular brackets denote classical traces over counterion (c) and 
macroion (m) coordinates. The mapping proceeds by formally tracing over the counterion coordinates: 

Z = <cxp(-/?tf cff )) m , (6) 

where H c ff — H m + F c is the effective one-component Hamiltonian and 

F c = -k B Tln (cxp [-0(H e + H mc )]) c (7) 

is the free energy of a nonuniform gas of counterions in the presence of the fixed macroions. Within pertur- 
bation theory [fja . l69l | , the counterion free energy can be expressed as 

F C = F + J dX =Fo + J dA {Hmc)x ' (8) 
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where F C (X) = -fc s Tln (exp [-/3(H C + XH mc )]) c , F Q = F c (0) = ~k B T In (exp(-/3H c )) is the unperturbed 
counterion free energy in the case of uncharged (yet volume-excluding) macroions, ( } A denotes a counterion 
trace with the macroions charged to a fraction A of their full charge, and the A-integral charges up the 
macroions. After formally adding and subtracting the energy of a uniform compensating negative background 
E b , Eq. © becomes 

F c = F OC P + [ dX (H mc ) x - E b , (9) 
Jo 

where Focp = Fo+Ef, is the free energy of a classical one-component plasma (OCP) in the presence of neutral 
hard spheres. The background and counterions alike are excluded from the hard cores of the macroions and 
therefore occupy a free volume V = V(l — 77), where r\ = (Tr/6)(N m /V)a 3 is the macroion volume fraction. 



2. Response Theory 

To make practical use of the one-component mapping, the counterion free energy must be approximated, 
for which purpose response theory provides a powerful framework. Because it proves more convenient to 
manipulate Fourier components of densities and pair potentials, we first note that the macroion Hamiltonian 
[Eq. @] and macroion-counterion interaction [Eq. Q}] can be equivalcntly expressed as 

H m = # HS + — ^v mm {k)[p m {k)p m {-k) - N m ] (10) 



2V> 

k 



and 



H mc = x77 ^ ^mc(fc)pm(k);Qc(-k), (11) 
k 

where the Fourier transform and its inverse are defined as 

p m (k) = J dv Pm {v)e-^ (12) 

k 

Equation (|llf) makes evident that H mc depends, through p c (k), on the response of the counterion density 
to the macroion charge density. The counterion response can be approximated by expanding the ensemble- 
averaged induced counterion density in a functional Taylor series in powers of the dimcnsionlcss macroion 
potential |6^, 0, E3 > u ( r ) = ~P J dr' u mc (|r— r'|)p m (r'). Expanding about zero macroion charge (u = 0), 
the counterion density can be expressed, in Fourier space, as 



(/5 c (k)} = x( fc )«mc(fc)/O m (k) + -^j ^x'(k',k - k')v mc (k')v mc (\k - k'|) 

k' 

x p m (k')/3 m (k-k') + ---, fc^O, (14) 



where \ and x' are i respectively, the linear and first nonlinear response functions of the uniform OCP. The 
response functions are directly related to the structure of the OCP according to x(fc) = —(3n c S(k) and 
X'(k',k - k') = (/?V/2)S" (3) (k',k - k'), where 

SW(ki,...,k„_i) = -^(p c (k 1 )---p c (k n _ 1 )p c (-k 1 ...-k„_ 1 )) (15) 

is the OCP n-particle static structure factor 69], S(k) = S^(k), and n c = N c /V is the average density of 
counterions in the free volume. The first term on the right side of Eq. I|14|) . which is linear in v mc (k) and 
/5 m (k), represents the linear response approximation, while the higher-order terms are nonlinear corrections. 
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Combining Eqs. (JSJ, (fTT|l . and lfH|) . specifying the background energy as Eb = linH;^o{^^c"-c^cc(fc)/2}, 
isolating the fc = terms, and integrating over A, produces the counterion free energy to third order in the 
macroion density: 



F c = Fqcp + n c lim 



N m v mc (k) + ^v cc (k) 



-^J E X< ^ ^rncik)] 2 /3 m (k)/3 m (-k) 



+ ^77^ E E x'(k', -k - k')v mc (k)v mc (k')v mc (\k + k'|)/5 m (k)/5 m (k')p m (-k - k'). 



(16) 



k^O k' 



The terms in F c that are quadratic and cubic in p m (k) generate effective pair and triplet interactions, 
respectively, in the effective Hamiltonian. To demonstrate this, we first identify 



= X(k)[Vrnc(k)] 2 

as an effective pair interaction, induced by linear response of counterions |57ll58ll68j |. and 

r>$(k, k') = 2 X '(k', -k - k')v mc (k)v mc (k')v mc (\k + k'|) 

as an effective three-body interaction, induced by nonlinear counterion response. Combining Eqs. (|10|l and 
(|f 6L the effective Hamiltonian now can be recast in the form 



(17) 



(18) 



1 N m N m 



(19) 



where v^g (r) = v mm {r) + v^ d (r) and v^g (r, r') are the effective macroion pair and triplet potentials, re- 
spectively, and £ is a one-body volume energy, composed of all terms in H e s independent of macroion 
coordinates. The volume energy accounts for the counterion entropy and macroion-counterion interaction 
energy and contributes density-dependent terms to the total free energy that can influence thermodynamic 
properties, as discussed below (Sec. UVf l. 

Explicit expressions for the effective interactions are obtained by invoking the identities 



N m 2 

E v indd R - R ,D = E ^ d (k)p,n(k)p m (-k) + Hm <>£>(*) - NmvW(0) 



(20) 



and 



N„ 



k#0 



E «^(Ri-B,,Ri-H* 



^2 E E €s (k, k')[p m (k)p m (k')p m (-k k>) 

k k' 

- 3p m (k)p m (-k) + 2N m ]. 
The volume energy, E = Eq + AE, is the sum of the linear response approximation [HJ l5Sj 



'm„(2) 

2 

and the first nonlinear correction [5 



£0 = Focp + -yX^O) + ^™« c Hm 



QV' 2 



E^i^ko-^E^c^o) 



k,k' 



(21) 



(22) 



(23) 



(2) (2) (2) 

Similarly, the effective pair interaction, ir ff (r) = Uq (r) + Air ff (r) is the sum of the linear response ap- 
proximation [57], , vjf\r) = v mm (r) + v^ d (r), and the first nonlinear correction, Av'j (r), whose Fourier 
transform is 



A^ ) W = ^E^ ) ( k I k ')-^ , (k,0). 



(24) 



G 



It is important to note that nonlinear counterion response generates not only effective many-body interac- 
tions, but also corrections to the effective pair and one-body interactions. It is these corrections [Eqs. (|23[l 
and i|24[l] whose impact on phase behavior we examine below in Sec. IIVI Note that the final terms on the 
right sides of Eqs. I|22l) - (|24l) originate from the charge neutrality condition, which required special treatment 
of the k = terms in Eqs. I|16l) and (|20|l . A simple physical interpretation of microion response and its 
connection to microion-induced effective interactions between macroions is discussed in ref. l59l. 



3. Random Phase Approximation 



Further progress requires specifying the OCP response functions. To this end, we note first that the 
counterions are usually characterized by relatively small electrostatic coupling parameters, T = \ B /a c I, 
where X B = /3z 2 e 2 /e is the Bjerrum length and a c — (3/47rn c ) 1 / 3 is the counterion-sphere radius. In 
such weakly-coupled plasmas, short-range correlations are often weak enough to justify a random phase 
approximation (RPA) whereby the two-particle direct correlation function (DCF) is approximated 

by its exact asymptotic limit: c^ 2 \r) = — (3v cc (r) or <^ 2 \k) — —Aitf3z 2 e 2 /ek 2 . The OCP linear and first 
nonlinear response functions then take the simple analytical forms 

-(3n c 



X(k) 



and 



X'(k,k' 



k B T 
2n 2 . 



l+K 2 /k 2 



X(fc) x (fc')x(|k + k'|), 



(25) 



(26) 



where k = -J 4irn c z 2 e 2 / ek B T is the Debye screening constant (inverse screening length). Higher-order non- 
linear response leads to higher-order terms in the effective Hamiltonian [Eq. (|19f) ]. which are here neglected. 

Practical expressions for the effective interactions follow from specifying the macroion-counterion interac- 
tion inside the macroion core so as to minimize counterion penetration - a strategy similar to that of the 
pseudopotential theory of simple metals [TlL W% . The choice 



s (r) = - 



Zze k 
e(l + kcl) 



r < a, 



(27) 



ensures zero counterion penetration (p c (r) — 0, r < a) at the level of linear response [4!j, |57J, |58| and virtually 
eliminates counterion penetration in the case of nonlinear response |59| . Substituting the Fourier transform 
ofEq. 



AirZze 2 



<1 



)k 2 



cos(fca) + — sm(ka) 



(28) 



into Eqs. (fT7f>. ((TH1) . and l(^ |l -Pl )l then yields the effective interactions. 

Upon reintroducing salt ions as a second species of microion |5^ | , analytical expressions are obtained |59j | 
for the volume energy and the effective pair potential. The volume energy is the sum of the linear response 
approximation 



2„2 



E n = F, 



plasma 



N r . 



Z z e 



2e 1 



k B T {N+ - N_) 2 
2 N+ + N_ 



and the first nonlinear correction 

N m k B T (n + - n_ 



AE = 



G 



Z 2 K 3 n u 


f 1 1 


2 Z 3 K 6 


( e Ka \ 


8» 


V 1 + Ka ) 


(An) 2 


V 1 + na J 



Ei(3ko) 



(29) 



(30) 



where F p i asma = k B T[N + ln(n + A 3 ) + A_ ln(n_A 3 )] is the ideal-gas free energy of the plasma, n± = N±/V' 
and n M = N^/V' = n + + n_ = n c + 2n s are the microion number densities in the free volume, A is the 
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thermal wavelength of the microions, and 

'4nz 2 e 2 nA 1/2 _ / We 2 (N C + 2N S )^ 1/2 



ek B T J V ^bT V(l - Tj) 



(31) 



is the Debye screening constant, which depends on the total density of microions, adjusted for macroion 
excluded volume. The effective pair potential is the sum of 

(2 ), , Z 2 e 2 ( e Ka \ 2 e~ Kr 



<^ = — {TT^a) — (32) 
which is identical to the DLVO potential with a density-dependent screening constant, and 

Av$ (r) = h (r) — + h (r) — + h M — , r > a, (33) 
r r r 



where 



/l(r) = d [n{r - a) + 1 - e,- Ka } + C 2 [Ej (/s(r - a)) + Ei (3/ea) - Ex («a)] , (34) 

/2(r) = -C 2 Ei(3/s(r + a)), (35) 

/ 3 (r)=C 2 [E 1 (2/s(r + a))-Ei(2K(r-a))] ) (36) 

a '<=-ih)5!f!(^y, (37) 

D n M e \ 1 + rea / 

1 fn+ - n-) Z 3 e 2 n 3 ( e Ka x 3 



C2 = jr- v " 2 ; 71— ' ( 38 ) 

07r n.„ ze \ 1 + Ka / 



and 



Ei(x) = / du , jc>0, (39) 



is the exponential integral function. The effective three-body interaction can be computed from the gener- 
alizations of Eqs. (HBJ), IJ2SJ, and (J2HJ, with the result [5Q| 

«^(n - r 2 ,n - r 3 ) = -k B T (n+ ~ n ~ } / drpxdn-rDpxd^-rDpxdrs-rl), (40) 



where 

{Z k 2 e Ka e~ KT 

~4¥ I + Ka ~r~' r>a ' (41) 

0, r < a, 

is the density of counterions around an isolated macroion. 

It is important to establish the accuracy of the effective interactions predicted by the nonlinear response 
theory described above. In a direct comparison with ab initio simulations |76| . first-order nonlinear correc- 
tions were shown to quantitatively match effective pair energies |59l |. Nevertheless, the effective interactions 
predicted by response theory should be tested further, perhaps by comparisons with nonlinear Poisson- 
Boltzmann theory. 
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B. Thermodynamic Phase Behavior 

1. Variational Theory 

The effective electrostatic interactions predicted by response theory provide the basic input required by 
statistical mechanical theories and computer simulations of the effective one-component model of charged 
colloids. The one-component model is considerably simpler than the (multi-component) primitive model, 
and thus a practical alternative for investigating thermodynamic phase behavior and other bulk properties 
of many-particle systems. Here we input effective interparticle interactions into an approximate variational 
theory for the free energy. The Helmholtz free energy F separates naturally into three contributions: 

F(T, V, N m ,N s ) = F id {T, V, N m ) + F C *{T, V, N m , N s ) + E{T, V, N m , N s ), (42) 

where F^ = A r m fcsT[ln(n m A^ [ ) — 1] is the exact ideal-gas free energy of a uniform fluid of macroions of 
thermal wavelength A m , F cx is the excess free energy, which depends on effective intermacroion interactions, 
and E is the one-body volume energy. Note that fox and E depend on the average densities of both macroions 
and salt ions. 

To approximate the excess free energy, we apply a variational approach based on first-order thermodynamic 
perturbation theory, as in ref. |5f)j |. Given a decomposition of the effective pair potential into reference and 
perturbation potentials, 

«$(r)=t&V)+»St(r), (43) 

an upper bound on the excess free energy density, / ex = F cx /V, is provided by the Gibbs-Bogoliubov 
inequality |6^| 

1 „2 f , „ ./ x .(2) 



/ox < /ref + ^"m / dr Srcf(*>pcrtM> ( 44 ) 

where / re f and g ro f(r) are the excess free energy density and radial distribution function, respectively, of the 
reference system. The short-range-repulsive form of the effective pair potential naturally suggests a hard- 

(2) 

sphere (HS) reference system. Thus, t^ef ( r ) = v BS\ r \ a), the pair potential between hard spheres of effective 

(2) (2) 

diameter d, and Vp CIt (r) — u; ff (r), r > d. The effective HS diameter provides a variational parameter with 
respect to which the right side of Eq. I(44|) can be minimized to impose a least upper bound on the excess 
free energy: 

) - mm |/ HS (n m ,n s ; (i) + 27rn^ y" dr r 2 g HS (r, n m ; d)v^ (r, n m , n s )\ . (45) 

Here fns( n m, n s ;d) and 5Hs( r i n m ; d) are, respectively, the excess free energy density and radial distribution 
function of the HS reference fluid, which we approximate by the essentially exact Carnahan-Starling and 
Verlet-Weis analytical expressions |69|. In practice, the exponential decay of t> cff (r) with r ensures rapid 
convergence of the perturbation integral in Eq. 1)45(1 . justifying the further approximation that c/hsM = 1 
for r > 5d. The accuracy of the variational theory in predicting the equation of state has been confirmed by 
independent comparisons with Monte Carlo simulation data |50|, |7j| ■ 

2. Grand Potential and Phase Coexistence 

For a system at fixed temperature, volume, and number of macroions, in osmotic equilibrium with a 
salt reservoir at fixed salt chemical potential fi s , the appropriate thermodynamic potential (minimized at 
equilibrium) is the semi-grand potential, 

n(T, V, N m , p.) = F(T, V, N m , N s ) - fi s N s = -pV + ^ m N m , (46) 

where p is the bulk pressure and /j, m is the chemical potential of the pseudomacroions. More precisely [i m 
is the change in free energy - at constant T and V - upon adding a bare macroion and its Z/z neutralizing 
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counterions and /z s is the change in free energy upon adding a charge-neutral pair of salt ions. The semi-grand 
potential density is then given by 

u{T, n m , n s ) = fl/V = f(T, n m , n s ) - /j, s n s = -p + p. m n m , (47) 

where / = F/V is the total free energy density and n s = N s /V is the number density of salt ion pairs in 
the system. At constant T, the differential relation 

dfi(T, V, N m , n s ) = -pdV + n m dN m - N s dfi s , (48) 

yields the pressure 

—) =n (—\ -u, (49) 

9F/m„ „ ™ \ dn m ) „ ., 



and the macroion chemical potential 



dCl\ ( duo 



"— 'aid™ - sri • (50) 



p(l) 


= p( 2 ) 




= 




= /4 2) 



Equilibrium coexistence of bulk phases requires equality of pressure and of chemical potentials (of 
macroions and salt ions) in the two phases (1 and 2): 

(51) 
(52) 

4 r) , (53) 

where the superscript (r) denotes a reservoir quantity. Equality of pressure is equivalent to equality of 
osmotic pressure, H = p — p^ r \ i.e., the difference between the system and reservoir pressures. The osmotic 
pressure - a manifestation of the Donnan effect jjj] - vanishes in the dilute limit of zero colloid concentration. 

The coexistence conditions have simple geometrical interpretations. Equations H47(l - (|53(l describe a com- 
mon tangent, of slope [i m and intercept — p, to the curve of w{n m , /j. s ) vs. n m (constant fi s ), or equivalently 
a Maxwell equal-area construction. Specifically, the relations 



2 



du) = I d» m ^ m (n ro , ft) = ^'(n^ - n^') (54) 



and 

r2 r2 



d(n/N m ) = - / dv mP (v m ^ s ) = - P W{vW - v£>), (55) 



with v m = V/N m = l/n m , imply that constant-/i s curves of /i m (n m , ji s ) vs. n m and of p(v m ,fi s ) vs. v m 
enclose equal areas above and below the horizontal lines \x m = \i~m = Mm* and p = — p^\ respectively. 
Changes of curvature sufficient to allow common-tangent constructions on the semi-grand potential, and 
equal-area constructions on the chemical potential and pressure, imply phase coexistence. 

At low salt concentrations, the salt reservoir behaves as an ideal gas of ions, whose pressure and chemical 
potential are well approximated by 



P 



M = 2n<fh B T (56) 



and 



/4 r) = 2/c B 7Tn(nWA 3 ), (57) 



where ni^ is the reservoir number density of pairs of salt ions of thermal wavelength A. Note that A and 
A m are arbitrary, as they contribute to the semi-grand potential only terms that are linear in density, which 
do not affect the coexisting densities. 
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The phase diagram is computed as follows. For a given macroion density n m and salt chemical potential 
(i.e., reservoir salt density ni ), the system salt density n s is numerically determined [from Eq. (|53Jl ] as the 
solution of 

„ s = ( df{n d ™l ns) ) T = 2k B THn^A% (58) 

where / is the total free energy density, the excess part of which is given by Eq. I|45|l . In the case of linear 
response, Eq. I|58|l can be expressed in a somewhat more practical form by separating out and analytically 
evaluating the dominant volume energy contribution. Substituting Eqs. (|2~9")l and into Eq. ifS^jl then 
yields 

.31 , , , A 3n ZkX b n c i ( n c \ 2 t a ( df ex (n m ,n s 



fo, = H(n c + n s )Al + Hn s A") 9n T , 2 ~ + ~ + H 7 • (59) 

2(1 + Kay \riftj \ dn s J Hm 

The pressure and macroion chemical potential are next computed from Eqs. I|49(l and (|50() . Finally, the 
macroion and salt densities are varied to satisfy the remaining coexistence conditions [Eqs. I|51|l and (|52l) ]. 



IV. RESULTS AND DISCUSSION 



To investigate the influence of nonlinear microion screening on the phase behavior of deionized charged 
colloids, the variational theory (Sec. 1111 fj|l is used to compute the semi-grand potential, taking as input the 
effective interactions predicted by response theory (Sec. IIII Ajl . By performing a coexistence analysis and 
comparing the phase diagrams that result from linear and first-order nonlinear interactions, leading-order 
nonlinear effects are quantified. For simplicity, effective three-body interactions are here neglected, since 
these are always attractive 59] and thus would only promote phase separation. In this way, we isolate the 
main nonlinear corrections to the volume energy and effective pair potential and assess their impact on phase 
behavior. 

Numerical results are presented for the case of room-temperature aqueous suspensions (As = 0.72 nm) 
and monovalent counterions (z — 1). For several choices of macroion radius a, the effective macroion 
valence Z is set near the threshold for charge renormalization [73], Z ~ O(10)(a/ Figure |2 illustrates 
the effective pair and triplet potentials vs. macroion separation, with linear and nonlinear screening, for 
various sets of system parameters. The particular case of (a = 266 nm, Z = 1217) is included to permit 
direct comparison with ref. jsjj. While nonlinear screening generally softens repulsive pair interactions, 
the correction is relatively minor for the selected macroion diameters and valences. The effective triplet 
potential, shown for an equilateral triangle arrangement of three macroions, is always attractive and decays 
rapidly with increasing separation. In passing, we note that the triplet interactions that arise within response 
theory | 59| differ in definition from their counterparts in Poisson-Boltzmann theory [74L l75| . 

Figures 13 and 0] present predictions for the osmotic pressure II (equation of state) vs. volume fraction 
r\ at fixed reservoir salt concentration (or salt chemical potential /i s ). The variation of II with r/ 
is a diagnostic of thermodynamic stability, a negative slope signaling instability toward phase separation 
(see below). Figure |3] illustrates that, within the linearized theory, the system becomes unstable below a 
certain critical salt concentration. Figure0]demonstrates the sensitivity of the osmotic pressure to nonlinear 
screening, which originates mainly from the nonlinear correction to the volume energy. 

Figures and present the corresponding system salt concentration c s (in ^mol/liter) vs. volume fraction 
(at fixed /z s ). The monotonic decrease of c s with increasing rj follows from Eq. (|59|) and stems from an 
interplay between salt entropy and salt-macroion interactions. Entropy and excluded-volume interactions 

(r) . (r) 

alone would give a simple linear decline, c s = (1 — ryjc s , with a slope of — c s . However, salt-macroion 
electrostatic interactions tend to expel salt from the system, steepening the decline, while maintaining an 
approximate linear dependence over a considerable range of r/. As illustrated in Fig. nonlinear screening, 
which modifies the state dependence of the effective interactions, tends to lower the system salt concentration. 

Figure typifies the monotonic decrease of the effective hard-sphere diameter d, and increase of the Debye 
screening constant k, with increasing volume fraction at fixed /z s . Nonlinear screening evidently reduces both 
d and k. For the chosen parameters, the reduction appears modest, but is significant, given the sensitivity 
of the free energy to these parameters. 
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Figures |3J and El illustrate that, for sufficiently high macroion valence and low salt concentration, 
van der Waals loops emerge in the equation of state at fixed [i s - a direct signature of phase instability. 
The maximimum and minimum in the curve of osmotic pressure vs. volume fraction mark the vapor and 
liquid spinodal densities, respectively, between which the compressibility is negative and the uniform fluid 
is unstable with respect to phase separation [Fig. Et a )]- Correspondingly, an equal-area construction on 
the curve of osmotic pressure vs. inverse volume fraction [Fig, ffib)]. or of chemical potential vs. volume 
fraction [Fig. Etc)], yields the densities of the coexisting vapor and liquid phases. A scan over reservoir salt 
concentration (salt chemical potential) traces out the spinodal and binodal (coexistence) curves in the phase 
diagram. 

Figure ED presents the resulting fluid phase diagrams for highly deionized suspensions as predicted by 
variational theory with both linear and nonlinear effective interactions as input. In each case, above a 
critical salt concentration, the uniform fluid is thermodynamically stable. Below the critical point, the fluid 
separates into macroion-rich (liquid) and macroion-poor (vapor) bulk phases, the salt concentration playing 
a role analogous to temperature in the liquid-vapor separation of a simple one-component fluid. For the 
parameter regime investigated here, the density of the liquid phase is found to be always well below the 
threshold for freezing, estimated from the hard-sphere freezing criterion, rj(d/a) 3 ~ 0.49, with the charged 
colloids approximated as neutral hard spheres of effective diameter d. 

The tie lines in the phase diagrams of Fig. join corresponding points on the liquid and vapor binodals 
(and spinodals) and, if extended, intersect the r\ = axis at the respective reservoir salt concentrations. 
The fact that the tie lines all have essentially the same slope, independent of reservoir salt concentration, 
is a physical consequence of strong salt-macroion electrostatic interactions, as described by Eq. I|59|l . The 
influence of nonlinear response on the tie-line slopes is negligible for the parameters here investigated. 

The predicted phase separation of charged colloids is remarkable, considering that simple one-component 
systems, interacting via purely repulsive pair potentials, exhibit only a single fluid phase. Within the 
present theoretical framework, phase instability at low salt concentrations is driven by the strong density 
dependence of the effective interactions, chiefly the one-body volume energy in deionized suspensions. It 
should be emphasized that because the colloid and salt concentrations vary between the two phases, the 
density-dependent effective interactions also differ in the two phases. 

The unusual phase separation can be understood, more fundamentally, as the result of a classic compe- 
tition between entropy and energy. On one side of the balance, favoring a stable uniform fluid, are the 
configurational entropies of all ions, represented by the ideal-gas terms in Eqs. l|29"|) and JI2J|, and the pos- 
itive potential energy of macroion pair repulsion. On the other side is the (density-dependent) negative 
potential energy of macroion-counterion attraction [second term on the right side of Eq. 1|29[) ] , which favors 
a concentrated phase with counterions localized around, and thus strongly attracted to, the macroions. 

Within the "entropy vs. energy" view, the sensitivity of phase behavior to salt concentration becomes 
clearer. At salt concentrations low enough that screening is counterion-dominated and screening lengths are 
relatively long, the counterion distribution is so diffuse that counterion-macroion attraction is too weak to 
drive macroion aggregation. With increasing salt concentration, the screening length shortens, the counteri- 
ons become more localized around the macroions, and counterion-macroion attraction may - for sufficiently 
high macroion valence - overcome configurational entropy and macroion pair repulsion to drive phase sepa- 
ration. The resulting concentrated phase is energetically favored, the counterions being closer on average to 
the macroions, but entropically disfavored, since the microions (excluded by macroion cores) must occupy a 
smaller free volume. On the other hand, the dilute phase is energetically disfavored, the counterions tending 
to roam farther from the macroions, but is entropically favored, since the microions can explore a larger free 
volume. At salt concentrations high enough that screening is salt-dominated, the salt-ion entropy overwhelms 
the counterion-macroion interaction energy in the free energy and prevents macroion aggregation. 

Thermodynamic phase behavior qualitatively similar to that depicted in Fig. |U| has been predicted be- 
fore [HoL l54|. Compared with the results of van Roij et al. jH^I, based on essentially the same variational 
theory for free energies, but a linearized density-functional theory for effective interactions, the present 
theory predicts a somewhat larger unstable area in the phase diagram. This quantitative discrepancy re- 
sults mainly from different treatments of excluded-volume effects in the two approaches. In particular, 
the excluded-volume correction to the screening constant in response theory [1/(1 — rf) factor in Eq. I|31|) ] 
enhances microion screening and promotes phase instability. 
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V. SUMMARY AND CONCLUSIONS 

In summary, we have investigated the controversial issue of phase separation in deionized charge-stabilized 
colloidal suspensions by inputting effective electrostatic interactions from response theory into free energies 
from a thermodynamic variational theory. By considering both linear and first-order nonlinear approxi- 
mations for the effective pair potential and one-body volume energy, we have systematically assessed the 
influence of nonlinear screening on phase behavior. A coexistence analysis results in osmotic pressures 
[Figs. 13 El and|H] and phase diagrams [Fig.E] that clearly exhibit thermodynamic instability towards phase 
separation for sufficiently high macroion effective valences and low salt concentrations. 

For macroion sizes and effective valences within limits established by charge renormalization considerations, 
first-order nonlinear corrections to the effective interactions are relatively weak and can either enhance or 
diminish stability of the uniform fluid phase, depending on system parameters. In general, the higher 
the macroion surface charge density, the higher the critical salt concentration and the larger the area of 
the unstable region in the phase diagram. Our main conclusion is that, within the present model, nonlinear 
screening appears not to suppress pha se separation of deionized suspensions, contradicting conclusions drawn 
from previous studies |3 |M H H3 and raising hope that a similar phenomenon may yet be observed in 
simulations of the primitive model. 

In closing, three key approximations of the present approach deserve to be highlighted for further scrutiny. 
First, the neglect of higher-order nonlinear corrections to the effective interactions presumes that nonlinear 
effects are strongest at the one- and two-body levels. The finding that first-order nonlinear corrections do 
not qualitatively alter fluid phase behavior suggests that higher-order corrections are unlikely to have drastic 
consequences - for example, suppression of phase separation. Furthermore, the presumption of weak many- 
body effective interactions is consistent with the dominance of the volume energy in effective one-component 
models of simple metals j^, |?3 H3, El E3i but should be further checked for charged colloids. Second, 
the mean-field approximation for the response functions of the microion plasma assumes weakly correlated 
microions. Although usually considered reasonable for monovalent microions, this assumption can and should 
be checked by more accurately modeling the structure of the microion plasma. Finally, the assumption of 
fixed macroion valence neglects the dependence of the effective valence on colloid and salt densities. This 
interesting issue of coupling between the effective macroion charge and phase behavior is being examined by 
means of charge renormalization theory and will be the subject of a future paper. 
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FIG. 1: Models of charge-stabilized colloidal suspensions: (a) Primitive model of charged hard-sphere macroions, of 
effective valence Z and diameter ct, and microions (counterions, salt ions) suspended in a dielectric continuum, (b) 
Effective one-component model of pseudomacroions governed by effective interactions. 
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FIG. 2: Effective pair potential (r) vs. center-to-center separation r for fixed colloid volume fraction r\ = 0.05 
and various combinations of macroion diameter a, effective valence Z, and system salt concentration c s : (a) a = 100 
ran, Z = 500, c s = 50 ^M; (b) a = 266 nm, Z = 1217, c s = 10 ^M; (c) a = 500 nm, Z = 2000, c s = 10 /uM. 
Solid (dashed) curves are predictions of nonlinear (linear) response theory. Insets show corresponding effective triplet 
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FIG. 3: Linear-screening predictions for osmotic pressure II (in reduced units) vs. colloid volume fraction r\ for same 
combinations of macroion diameter a and valence Z as in Fig. gland various fixed reservoir salt concentrations 



(a) a = 100 nm, Z 
Z = 2000 



500, 



100, 200, 400 fiM; (b) a = 266 nm, Z = 1217, c{ r) = 40, 80, 160 pM; (c) a = 500 nm 



10, 20,40 [M. 



20 




0.05 0.1 

Volume Fraction r| 




CO 



0.05 0.1 

Volume Fraction r| 



21 



GO 



CO 

c 

CD 

13 
CO 
CO 
CD 

CL 

O 
-i— » 
O 

E 

CO 

O 







1 1 1 1 1 

(c) g=500 nm 




Z=2000 


/ ■ 


c s (r) =20 |iM 




Linear 


/ : 


Nonlinear / 


/ 

/ ■ 

/ 




/ 


— * " ^ — * 


/ 

/ 

/ 

— 


i 









0.05 

Volume Fraction r| 



0.1 



FIG. 4: Osmotic pressure II (in reduced units) vs. colloid volume fraction 77 for same combinations of macroion 
diameter a and valence Z as in Fig. and fixed reservoir salt concentration ci r ^: (a) a — 100 nm, Z = 500, 
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curves are predictions of nonlinear (linear) response theory. 
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FIG. 7: Effective hard-sphere diameter d (units of macroion diameter a) and Debye screening constant k (inset) vs. 
colloid volume fraction rj, at fixed reservoir salt concentration c[ r \ for (a) a = 100 nm, Z — 500, ci r ' = 50 fjM; (b) 
a = 266 nm, Z = 1217, ci r) = 10 ^M; (c) o = 500 nm, Z = 2000, c{ r) = 10 £tM. Solid (dashed) curves are predictions 
of nonlinear (linear) response theory. 
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FIG. 8: Linear-screening prediction for (a) osmotic pressure II vs. colloid volume fraction r\, (b) II vs. I/77, and 
(c) colloid chemical potential /j, m (shifted by arbitrary constant) vs. 77 for macroion diameter a = 100 nm, valence 
Z = 500, and reservoir salt concentration ci r ' = 350 fjM. In panels (a) and (b), dotted vertical lines at maximum and 
minimum of II indicate spinodal densities at boundaries of unstable region. In panels (b) and (c), dashed vertical 
lines indicate coexisting densities on the fluid binodal, illustrating the Maxwell equal-area construction. 
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FIG. 9: Fluid phase diagrams for aqueous suspensions of charged colloids at room temperature (As = 0.72 nm) 
with monovalent microions and various macroion diameters and effective valences: (a) a = 100 nm, Z = 500; (b) 
a — 266 nm, Z = 1217; (c) a — 500 nm, Z — 2000. Solid (long-dashed) curves represent predictions for binodals from 
nonlinear (linear) response theory. Short-dashed curves represent predictions for spinodals (linear response only). 
Circular symbols denote critical points. Tie lines join corresponding points on liquid and vapor branches of binodals. 



